Causal viscous hydrodynamics in 2+1 dimensions for relativistic heavy-ion collisions 
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We explore the effects of shear viscosity on the hydrodynamic evolution and final hadron spectra of 
Cu+Cu collisions at ultrarelativistic collision energies, using the newly developed (2+l)-dimensional 
viscous hydrodynamic code VISH2+1. Based on the causal Israel-Stewart formalism, this code de- 
scribes the transverse evolution of longitudinally boost-invariant systems without azimuthal symme- 
try around the beam direction. Shear viscosity is shown to decelerate the longitudinal and accelerate 
the transverse hydrodynamic expansion. For fixed initial conditions, this leads to a longer quark- 
gluon plasma (QGP) lifetime, larger radial flow in the final state, and flatter transverse momentum 
spectra for the emitted hadrons compared to ideal fluid dynamic simulations. We find that the 
elliptic flow coefficient V2 is particularly sensitive to shear viscosity: even the lowest value allowed 
by the AdS/CFT conjecture r\j s > l/4n suppresses vi enough to have significant consequences for 
the phenomenology of heavy-ion collisions at the Relativistic Heavy Ion Collider. A comparison 
between our numerical results and earlier analytic estimates of viscous effects within a blast-wave 
model parametrization of the expanding fireball at freeze-out reveals that the full dynamical theory 
leads to much tighter constraints for the specific shear viscosity rj/s, thereby supporting the notion 
that the quark-gluon plasma created at RHIC exhibits almost "perfect fluidity". 

PACS numbers: 25.75.-q, 12.38.Mh, 25.75.Ld, 24.10.Nz 



I. INTRODUCTION 

Hydrodynamics is an efficient tool to describe the ex- 
pansion of the fireballs generated in relativistic heavy-ion 
collisions. As a macroscopic description that provides 
the 4-dimcnsional space-time evolution of the energy- 
momentum tensor T^ u (x) it is much less demanding than 
microscopic descriptions based on kinetic theory that 
evolve the (on-shell) distribution function f(x,p) in 7- 
dimcnsional phase-space. 

Ideal fluid dynamics is even more efficient since it re- 
duces the number of independent fields needed to de- 
scribe the symmetric energy-momentum tensor from 10 
to 5: the local energy density e(x), pressure p(x) and 
the normalized flow 4- velocity u^(x) (which has 3 inde- 
pendent components). The equation of state (EOS) p(e) 
provides a further constraint which closes the set of four 
equations d u T» v (x) = 0. 

Ideal fluid dynamics is based on the strong assump- 
tion that the fluid is in local thermal equilibrium and 
evolves isentropically. While local momentum isotropy 
in the comoving frame is sufficient for a unique decom- 
position of the energy-momentum tensor in terms of e, p 
and it M , it does not in general guarantee a unique rela- 
tionship p(e). Gcnerically, the equation of state p(e) (a 
key ingredient for closing the set of hydrodynamic equa- 
tions) becomes unique only after entropy maximization, 
i.e. after a locally thermalized state, with Maxwcllian 
(or Bosc-Einstein and Fcrmi-Dirac) momentum distri- 
butions in the comoving frame, has been reached. For 



* Correspond to song@mps.ohio-state.edu 



this assumption to be valid, the miscroscopic collision 
time scale must be much shorter than the macroscopic 
evolution time scale. Since the fireballs created in rel- 
ativistic heavy-ion collisions are small and expand very 
rapidly, applicability of the hydrodynamic approach has 
long been doubted. 

It came therefore as a surprise to many that the bulk 
of the matter produced in Au+Au collisions at the Rela- 
tivistic Heavy Ion Collider (RHIC) was found to behave 
like an almost ideal fluid. Specifically, ideal fluid dy- 
namic models correctly reproduce the hadron transverse 
momentum spectra in central and semi-peripheral colli- 
sions, including their anisotropy in non-central collisions 
given by the elliptic flow coefficient V2(pt) and its depen- 
dence on the hadron rest mass, for transverse momenta 
up to about 1.5-2 GeV/c [l[ which covers more than 99% 
of the emitted particles. This observation has led to the 
conclusion that the quark-gluon plasma (QGP) created 
in RHIC collisions thermalizes very fast and must there- 
fore be strongly (non-perturbatively) interacting Q , giv- 
ing rise to the notion that the QGP is a strongly coupled 
plasma [1, 0, HI that behaves like an almost perfect fluid 

At RHIC energies, the almost perfect ideal fluid 
dynamical description of experimental data gradually 
breaks down in more peripheral collisions, at high trans- 
verse momenta, and at forward and backward rapidities; 
at lower energies it lacks quantitative accuracy even in 
the most central collisions at midrapidity [3]. Most of 
these deviations from ideal fluid dynamical predictions 
can be understood as the result of strong viscous effects 
during the late hadronic stage of the fireball expansion 
[3] after the QGP has hadronized. As the initial energy 
density of the fireball decreases, the dissipativc dynam- 
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ics of the hadronic stage takes on increasing importance, 
concealing the perfect fluidity of any quark-gluon plasma 
possibly created at the beginning of the collision. How- 
ever, as also pointed out in [3], persisting uncertainties 
about the initial conditions in heavy-ion collisions leave 
room for a small amount of viscosity even during the 
early QGP stage. Furthermore, the observed deviations 
of the elliptic flow parameter v^ipr) at large px even 
in the largest collision systems at the highest available 
collision energies are consistent with viscous effects dur- 
ing the early epoch of the fireball [l5|, [l6[ . During this 
epoch, the matter is so dense and strongly interacting 
that a microscopic description based on classical kinetic 
theory of on-shell partons [l5[ may be questionable. We 
therefore develop here a dissipative generalization of the 
macroscopic hydrodynamic approach, viscous relativistic 
fluid dynamics. 

The need for such a framework is further highlighted 
by the recent insight that, due to quantum mechanical 
uncertainty [13], no classical fluid can have exactly van- 
ishing viscosity (as is assumed in the ideal hydrodynamic 
approach). Even in the limit of infinitely strong coupling, 
the QGP must hence maintain a nonzero viscosity. Re- 
cent calculations of the shear viscosity to entropy ratio 
(the "specific shear viscosity" rj/s) in a variety of con- 
formal gauge field theories which share some properties 
with QCD, using the AdS/CFT correspondence, suggest 
a lower limit of 2 = -L [ll, [H, H3] - This is much smaller 
than the value obtained from weak coupling calculations 
in QCD [2l[ (although close to a recent first numerical 
result from lattice QCD [13]) and more than an order of 
magnitude below the lowest measured values in standard 
fluids [IH . Some alternative ideas how small effective vis- 
cosities could be generated by anomalous effects (chaotic- 
ity) in anisotropically expanding plasmas [23| or by neg- 
ative eddy viscosity in 2-dimensional turbulent flows [24| 
have also been proposed. 

Initial attempts to formulate relativistic dissipative 
fluid dynamics as a relativistic generalization of the 
Navier-Stokes equation [25|, [26| ran into difficulties be- 
cause the resulting equations allowed for acausal signal 
propagation, and their solutions developed instabilities. 
These difficulties are avoided in the "2nd order formal- 
ism" developed 30 years ago by Israel and Stewart (27| 
which expands the entropy current to 2nd order in the 
dissipative flows and replaces the instantaneous identifi- 
cation of the dissipative flows with their driving forces 
multiplied by some transport coefficient (as is done in 
Navier-Stokes theory) by a kinetic equation that evolves 
the dissipative flows rapidly but smoothly towards their 
Navier-Stokes limit. This procedure eliminates causality 
and stability problems at the expense of numerical com- 
plexity: the dissipative flows become independent dy- 
namical variables whose kinetic equations of motion are 
coupled and must be solved simultaneously with the hy- 
drodynamic evolution equations. This leads effectively to 
more than a doubling of the number of coupled partial 
differential equations to be solved [28j . 



Only recently computers have become powerful enough 
to allow efficient solution of the Israel-Stewart equations. 
The last 5 years have seen the development of numeri- 
cal codes which solve these equations (or slight variations 
thereof [13, H, H, HO, EH ) numerically, for systems with 
boost-invariant l ong itudinal expansion and transverse ex- 
pansion in zero l29l. l3fl ], one [3(| H3, l33l |34| and two di- 
mensions [35L l3€3l. |3TL138 I (see also Ref. [39j] for a numerical 
study of the relativistic Navier-Stokes equation in 2+1 
dimensions). The process of verification and validation 
of these numerical codes is still ongoing: While different 
initial conditions and evolution parameters used by the 
different groups of authors render a direct comparison of 
their results difficult, it seems unlikely that accounting 
for these differences will bring all the presently available 
numerical results in line with each other. 

We here present results obtained with an indepen- 
dently developed (2+l)-dimcnsional causal viscous hy- 
drodynamic code, VISH2+1 |4(j. While a short ac- 
count of some of our main findings has already been 
published [13], we here report many more details, in- 
cluding extensive tests of the numerical accuracy of the 
code: We checked that (i) in the limit of vanishing vis- 
cosity, it accurately reproduces results obtained with the 
(2+l)-d ideal fluid code AZHYDRO [H; (ii) for homo- 
geneous transverse density distributions (i.e. in the ab- 
sence of transverse density gradients and transverse flow) 
and vanishing relaxation time it accurately reproduces 
the known analytic solution of the relativistic Navier- 
Stokes equation for boost-invariant 1-dimensional lon- 
gitudinal expansion [42| : (iii) for very short kinetic re- 
laxation times our Israel-Stewart code accurately repro- 
duces results from a separately coded (2+l)-d relativis- 
tic Navier-Stokes code, under restrictive conditions where 
the latter produces numerically stable solutions; and (iv) 
for simple analytically parametrized anisotropic velocity 
profiles the numerical code correctly computes the ve- 
locity shear tensor that drives the viscous hydrodynamic 
effects. 

In its present early state, and given the existing open 
questions about the mutual compatibility of various nu- 
merical results reported in the recent literature, we be- 
lieve that it is premature to attempt a detailed compar- 
ison of VISH2+1 with experimental data, in order to 
empirically constrain the specific shear viscosity of the 
QGP. Instead, we concentrate in this paper on describ- 
ing and trying to understand the robustness of a variety 
of fluid dynamical effects generated by shear viscosity in 
a relativistic QGP fluid. We report here only results for 
Cu+Cu collisions, with initial entropy densities exceed- 
ing significantly those that can be reached in such colli- 
sions at RHIC. The reasons for doing so arc purely tech- 
nical: Initially our numerical grid was not large enough 
to accomodate Au+Au collision fireballs with sufficient 
resolution, and while this restriction has been lifted in 
the meantime, a large body of instructive numerical re- 
sults had already been accumulated which would have 
been quite expensive to recreate for the Au+Au system. 
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The unrealistic choice of initial conditions was driven by 
the wish to allow for a sufficiently long lifetime of the 
QGP phase even in peripheral Cu+Cu collisions such 
that shear viscous effects on elliptic flow arc still dom- 
inated by the quark-gluon plasma stage. The main goals 
of the present paper are: (i) to quantitatively estab- 
lish shear viscous effects on the evolution of the energy 
and entropy density of the flow profile, source eccentric- 
ity, and total momentum anisotropy, on the final hadron 
spectra, and on the elliptic flow in central and non-central 
heavy-ion collisions, under the influence of different equa- 
tions of state; and (ii) to explore in detail and explain 
physically how these effects arise, trying to extract gen- 
eral rules and generic features which should also apply for 
other collision systems and collision energies. We note 
that recent calculations for Au+Au collisions [43| have 
shown that viscous effects are somewhat bigger in the 
smaller Cu+Cu studied here than in the larger Au+Au 
system for which the largest body of experimental data 
exists. The reader must therefore apply caution when 
trying to compare (in his or her mind's eye) our results 
with the well-known RH1C Au+Au data. 

The paper is organized as follows: Section [TT] gives a 
brief review of the Israel-Stewart formalism for causal 
relativistic hydrodynamics for dissipative fluids, lists the 
specific form of these equations for the (2+l)-dimensional 
evolution of non-central collision fireballs with boost- 
invariant longitudinal expansion which are solved by 
VISH2+1, and details the initial consditions and the 
equation of state (EOS) employed in our calculations. In 
Section IIHI we report results for central (b = 0) Cu+Cu 
collisions. Section IIVI gives results for non-central colli- 
sions, including a detailed analysis of the driving forces 
behind the strong shear viscous effects on elliptic flow 
observed by us. In Section [V] we explore the influence of 
different initializations and different relaxation times for 
the viscous shear pressure tensor on the hydrodynamic 
evolution and establish the limits of applicability for vis- 
cous hydrodynamics in the calculation of hadron spectra. 
Some technical details and the numerical tests performed 
to verify the accuracy of the computer code are discussed 
in Appendices lAlDl and in Appendix [E] we compare our 
hydrodynamic results with analytical estimates of shear 
viscous effects by Teaney [lj| that were based on Navier- 
Stokes theory and a blast-wave model parametrization of 
the fireball. 



II. ISRAEL-STEWART THEORY OF CAUSAL 
VISCOUS HYDRODYNAMICS 

In this section, we review briefly the 2nd order Israel- 
Stewart formalism for viscous relativistic hydrodynamics 
and the specific set of equations solved by VISH2+I for 
anisotropic transverse expansion in longitudinally boost- 
invariant fireballs. Details of the derivation can be found 
in [28|, w ith a small correction pointed out by Baier et 
al. in [31]. For simplicity, and in view of the intended ap- 



plication to RHIC collisions whose reaction fireballs have 
almost vanishing net baryon density, the discussion will 
be restricted to viscous effects, neglecting heat conduc- 
tion and working in the Landau velocity frame [2(| . 

A. Basics of Israel-Stewart theory 

The general hydrodynamic equations arise from the 
local conservation of energy and momentum, 

d^(x) = 0, (1) 

where the energy-momentum tensor is decomposed in the 
form 

T» v = en^u" - (p+n)A^ +tt^. (2) 

Here e and p are the local energy density and thermal 
equilibrium pressure, and is the (timelike and nor- 
malized, u^u^ = I) 4-velocity of the energy flow. II is 
the bulk viscous pressure; it combines with the thermal 
pressure p to the total bulk pressure. In Eq. © it is 
multiplied by the projector A M " = g^ u — u^u v transverse 
to the flow velocity, i.e. in the local fluid rest frame the 
bulk pressure is diagonal and purely spacelike, (p + U)5ij. 
n^" is the traceless shear viscous pressure tensor, also 
transverse to the 4-velocity (tt^u^ = 0) and thus purely 
spatial in the local fluid rest frame. 

For ideal fluids, II and ■n^ vanish, and the only dy- 
namical fields are e(x), p(x) and u^i^x), with e and p 
related by the equation of state p(e). In dissipative flu- 
ids without heat conduction, Id and the 5 independent 
components of w^" enter as additional dynamical vari- 
ables which require their own evolution equations. In 
relativistic Navier-Stokes theory, these evolution equa- 
tions degenerate to instantaneous constituent equations, 

n=-CVu, tt^ = 2i 1 V^u u) , (3) 

which express the dissipative flows II and 7r Ml/ directly 
in terms of their driving forces, the local expansion 
rate 6 = V-u and velocity shear tensor a^ v = V^u v \ 
multiplied by phenomenological transport coefficients 
Ci V > (the bulk and shear viscosity respectively). Here 
V = A^ v d v is the gradient in the local fluid rest frame, 
and V^u l/ > = ±(V^w l/ +V 1 '?^) - ±(V-w)A^, showing 
that, like tt^, the velocity shear tensor is traceless and 
transverse to u^. The instantantancous identification © 
leads to causality problems through instantaneous signal 
propagation, so that this straightforward relativistic gen- 
eralization of the Navier-Stokes formalism turns out not 
to be a viable relativistic theory. 

The Israel-Stewart approach [27| avoids these problems 
by replacing the instantaneous identifications (J3|) with 
the kinetic evolution equations 

du = — (n + cv-u), (4) 

m 

Dir^ = —(ir> lv -2rf7 {ll u v) ) 

-( U ^TT Va +U V TT^)Du a , (5) 
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where D = u^d^ is the time derivative in the local fluid 
rest frame, and the last term in the bottom equation 
ensures that the kinetic evolution preserves the trace- 
lessness and transversality of 7r M1/ [4j|. rn and are 
relaxation times and related to the 2nd order expansion 
coefficients in the entropy current [13, H^j ■ The fact that 
in the Israel-Stewart approach the dissipativc flows II 
and 7r M1/ no longer respond to the corresponding thermo- 
dynamic forces V • u and V^«"' instantaneously, but on 
finite albeit short kinetic time scales, restores causality 
of the theory [27| . 

We should caution that the form of the kinetic evolu- 
tion equations (|4I5|I is not generally agreed upon, due to 
unresolved ambiguities in their derivation [27l.[28l.[29l[30l. 
l3lL [321 . l35l . [H, HJ, HI] . We will here use the form given in 
Eqs. (|4I5[) and comment on differences with other work 
when we discuss our results. 

In the following calculations we further simplify the 
problem by neglecting bulk viscosity. Bulk viscosity van- 
ishes classically for a system of massless partons, and 
quantum corrections arising from the trace anomaly of 
the energy-momentum tensor are expected to be small, 
rendering bulk viscous effects much less important than 
those from shear viscosity. This expectation has been 
confirmed by recent lattice calculations 0, [4?J which 
yield very small bulk viscosity in the QGP phase. The 
same calculations show, however, a rapid rise of the 
bulk viscosity near the quark-hadron phase transition 
[13] , consistent with earlier predictions [H, [49[ . In the 
hadronic phase it is again expected to be small [48|. We 
leave the discussion of possible dynamical effects of bulk 
viscosity near the quark-hadron phase transition for a 
future study. Bulk viscous pressure effects can be easily 
restored by substituting p — > p+Tl everywhere below and 
adding the kinetic evolution equation (U) for n. 

B. Viscous hydrodynamics in 2+1 dimensions 

In the present paper we eliminate one of the three 
spatial dimensions by restricting our discussion to lon- 
gitudinally boost-invariant systems. These are conve- 
niently described in curvilinear x m = (r, x, y, 77) coordi- 
nates, where r = V t 2 —z 2 is the longitudinal proper time, 
^ = 5 m (f__§) is the space-time rapidity, and (x,y) are 
the usual Cartestsian coordinates in the plane transverse 
to the beam direction z. In this coordinate system, the 
transport equations for the full energy momentum tensor 
T^ v are written as [HI 

d T f TT + d x (v x f TT ) + dy(v y f TT ) = s TT , 

d T f TX + d x (v x T™) + d y (v y T™) = S™, (6) 

d T T Ty + d X (V x T Ty ) + d y (V y T TV ) = S^. 

Here f mn = T(T m ™+7r m "), T m ™ = eu m u n -pA mn be- 
ing the ideal fluid contribution, u m = (u T , u x , u v , 0) = 
7x(l) v x ,v y) 0) is the flow profile (with 71 = , 1 2 _ 2 ), 



and g mn = diag(l, —1,-1, — 1/r 2 ) is the metric tensor for 
our coordinate system. The source terms S mn on the 
right hand side of Eqs. © are given explicitly as 

S TT = -p-r 2 ^ -Td x {pv x +w XT -v x w TT ) 
- Tdy(pv y +ir yT -v y ir TT ) 

~ -p - tV''' - rd x (pv x ) - Tdy(pv y ), (7) 
S TX = - T d x (p +7r xx -v x ir TX )-Tdy(ir xy -v y ir TX ) 

« -Td x (p+« xx ), (8) 

S Ty = -Td x (TT Xy -V x TT Ty ) -Tdy( P +TT yy ~Vy7r Ty ) 

W -Td y ( P +TT yy ). (9) 

We will see later (see the right panel of Fig. [TJ] below) 
that the indicated approximations of these source terms 
isolate the dominant drivers of the evolution and provide 
a sufficiently accurate quantitative understanding of its 
dynamics. 

The transport equations for the shear viscous pressure 
tensor are 

(d T + v x d x + v y d y )n mn = —(n mn -2 V a mn ) 

7_l_T.Tr 

-(u m n n 3 +u n n™)(d T + v x d x + v y d y y. (10) 

The expressions for a mn and n rnn are found in 
Eqs. (jAllA2|) ; they differ from n mn in Eqs. (7EJJ and a mn 
given in Ref. [-281 ] by a Jacobian r 2 factor in the (7777)- 
component: 7r ?,I, =t 2 tt t ' 71 , a m =T 2 a r,n . This factor arises 
from the curved metric where the local time derivative 
D = u m d m must be evaluated using covariant derivatives 
d m [50| . Since u v = 0, no such extra Jacobian terms arise 
in the derivative Dv? in the second line of Eq. (fTO)) . 

The algorithm requires the propagation of 7r rr , ~k tx ', 
ir Ty , and ir vv with Eq. Jp even though one of the first 
three is redundant (see [33] and Appendix [B]) . In ad- 
dition, we have chosen to evolve several more, formally 
redundant components of 7r m ™ using Eq. (J3J), and to use 
them for testing the numerical accuracy of the code, by 
checking that the transversality conditions u m 7T mn = 
and the tracelessness 7r^ = are preserved over time. We 
find them to be satisfied with an accuracy of better than 
1 — 2% everywhere except for the fireball edge where the 
Tr mn are very small and the error on the transversality 
and tracelessness constraints can become as large as 5%. 



C. Initial conditions 

For the ideal part TJ" T , T™, of the energy momen- 
tum tensor we use the same initialization scheme as for 
ideal hydrodynamics. For simplicity and ease of com- 
parison with previous ideal fluid dynamical studies we 
here use a simple Glauber model initialization with zero 
initial transverse flow velocity where the initial energy 
density in the transverse plane is taken proportional to 
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the wounded nucleoli density [44 
eo{x, y, b) = Kn WN (x, y; b) 



the approach to Navier-Stokes theory, we also did a few 
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Here a is the total inelastic nucleon-nucleon cross 
section for which we take er = 40mb. Ta,b is the 
nuclear thickness function of the incoming nucleus 
A or B. T A (x,y) = ^ oa dzp A {x,y, z); p A (x,y,z) is 
the nuclear density given by a Woods-Saxon profile: 
Pa{t) = i +CKp [ (r"-R A )/z } ■ For Cu+Cu collisions we take 
i? Cu = 4.2 fm, £ = 0.596 fm, and p Q = 0.17fm- 3 . The pro- 
portionality constant K does not depend on collision ccn- 
trality but on collision energy; it fixes the overall scale 
of the initial energy density and, via the associated en- 
tropy, the final hadron multiplicity to which it must be 
fitted as a function of collision energy. We here fix it 
to give eo = e(0, 0; 6=0) = 30 GeV/fm 3 for the peak en- 
ergy density in central Cu+Cu collisions, at an initial 
time To for the hydrodynamic evolution that we set as 
To = 0.6fm/c. As already mentioned in the Introduction, 
this exceeds the value reached in Cu+Cu collisions at 
RHIC (it would be appropriate for central Au+Au col- 
lisions at v / * = 200AGeV It ensures, however, a 
sufficiently long lifetime of the QGP phase in Cu+Cu 
collisions that most of the final momentum anisotropy is 
generated during the QGP stage, thereby permitting a 
meaningful investigation of QGP viscosity on the elliptic 
flow generated in the collision. 

Lacking a microscopic dynamical theory for the early 
pre-equilibrium stage, initializing the viscous pressure 
tensor 7r m ™ requires some guess-work. The effect of differ- 
ent choices for the initial Tr mn on viscous entropy produc- 
tion during boost-invariant viscous hydrodynamic evolu- 
tion without transverse expansion was recently studied 
in [5l|. We will here explore two options: (i) we set 
7T™" = initially [36[ ; or (ii) we assume that at time To 
one has it™ 1 = 2r)<j™ n where the shear tensor a™™ is cal- 
culated from the initial velocity profile u m — (1, 0, 0, 0) 
[33l . l35l l38l [39| . The second option is the default choice 
for most of the results shown in this paper. It gives 
t 2 ttq^ = —2itq x = ~2TtQ V = i.e. a negative contri- 

bution to the longitudinal pressure and a positive contri- 
bution to the transverse pressure. 

We here present results only for one value of the spe- 
cific shear viscosity, | = j- ~ 0.08, corresponding to its 
conjectured lower limit [la ]. The kinetic relaxation time 
Ttt will be taken as t, = J except were otherwise men- 
tioned. This value is half the one estimated from classi- 
cal kinetic theory for a Boltzmann gas of non-interacting 
massless partons [27], Hl[ - we did not use the twice larger 
classical value because it led to uncomfortably large vis- 
cous pressure tensor components n mn at early times, 
caused by large excursions from the Navier-Stokes limit. 
To study the sensitivity to different relaxation times and 
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in Section [VB] 



D. EOS 

Through its dependence on the Equation of State 
(EOS), hydrodynamic flow constitutes an important 
probe for the existence and properties of the quark- 
hadron phase transition which softens the EOS near T c . 
To isolate effects induced by this phase transition from 
generic features of viscous fluid dynamics we have per- 
formed calculations with two different equations of state, 
EOS I and SM-EOS Q. They are described in this sub- 
section. 

EOS I models a non-interacting gas of massless quarks 
and gluons, with p=^e. It has no phase transition. 
Where needed, the temperature is extracted from the 
energy density via the relation e= (16 + n^-Nf) fo^icF' 
corresponding to a chemically equilibrated QGP with 
Nf = 2.5 effective massless quark flavors. 




e(GeV/fm ) 

FIG. 1: The equations of state EOS Q (dashed line) and SM- 
EOS Q ("modified EOS Q", solid line). 

SM-EOS Q is a smoothed version of EOS Q [3 
which connects a noninteracting QGP through a first or- 
der phase transition to a chemically equilibrated hadron 
resonance gas. In the QGP phase it is defined by the 
relation p=|e— |S (i.e. c s = |f = |)- The vacuum 
energy (bag constant) B 1 ' 4 = 230 MeV is a parameter 
that is adjusted to yield a critical temperature T c = 
164 MeV. The hadron resonance gas below T c can be 
approximately characterized by the relation p = 0.15e 
(i.e. ci? = 0.15) fUj]. The two sides are matched through 
a Maxwell construction, yielding a relatively large la- 
tent heat Aei a t = 1.15 GeV/fm 3 . For energy densities 
between en = 0.45 GeV/fm 3 and eQ = 1.6 GeV/fm 3 one 
has a mixed phase with constant pressure (i.e. c^=0). 
The discontinuous jumps of c 2 s from a value of 1/3 to 



6 



at eqj and back from to 0.15 at en generate prop- 
agating numerical errors in VISH2+1 which grow with 
time and cause problems. We avoid these by smooth- 
ing the function c 2 (e) with a Fermi distribution of width 
Se = 0.1 GeV/fm 3 centered at e = eQ and another one of 
width 6e = 0.02 GeV/fm 3 centered at e = en- Both the 
original EOS Q and our smoothed version SM-EOS Q 
are shown in Figure [TJ A comparison of simulations us- 
ing ideal hydrodynamics with EOS Q and SM-EOS Q is 
given in Appendix ID II It gives an idea of the magnitude 
of smoothing effects on the ideal fluid evolution of elliptic 
flow. 

Another similarly smoothed EOS that matches a 
hadron resonance gas below T c with lattice QCD data 
above T c has also been constructed. Results using this 
lattice based EOS will be reported elsewhere. 



E. Freeze-out procedure: Particle spectra and V2 

Final particle spectra are computed from the hydro- 
dynamic output via a modified Cooper-Frye procedure 
[52j . We here compute spectra only for directly emitted 
particles and do not include feeddown from resonance 
decay after freeze-out. We first determine the freeze-out 
surface £(x), by postulating (as common in hydrody- 
namic studies) that freeze-out from a thermalized fluid 
to free-streaming, non-interacting particles happens sud- 
denly when the temperature drops below a critical value. 
As in the ideal fluid case with EOS Q [44j we choose 
Tdcc = 130 McV. The particle spectrum is then computed 
as an integral over this surface, 



E- 



d 3 p (2tt) 3 



p- d 3 a(x) fi(x,p) 



(12) 



(2tt) e 



p ■ d 3 a(x) [feq,i(x,p) + 8f l (x,p)} 



where gt is the degeneracy factor for particle species i, 
d 3 a fl (x) is the outward-pointing surface normal vector on 
the decoupling surface E(x) at point x, 



p ■ d a(x) = Utit cosh(y— rf) 
x Tf(r)rdrd<pdrj 



(13) 



(with r = (x, y) = (r cos cf>, r sin (j>) denoting the transverse 
position vector), and fi(x,p) is the local distribution 
function for particle species i, computed from hydro- 
dynamic output. Equation (|12|) generalizes the usual 
Cooper-Frye prescription for ideal fluid dynamics [52| by 
accounting for the fact that in a viscous fluid the local 
distribution function is never exactly in local equilibrium, 
but deviates from local equilibrium form by small terms 
proportional to the non-equilibrium viscous flows (l6ll3ll|. 
Both contributions can be extracted from hydrodynamic 
output along the freeze-out surface. The equilibrium con- 
tribution is 



feq,i{Pi X ) = /< 



cq.z 



p-u(x) 



1 



e p-u(x)/T(x) ± I ' 



(14) 



where the exponent is computed from the tem- 
perature T(x) and hydrodynamic flow velocity 
= 7^ (cosh?/, v± cos </>„, v± sin</>„, sinhr?) along the 
surface £(x): 

p ■ u(x) = j±[rriT cosh(y— 77) — prv± cos(4> p —(j) v )]. (15) 



Here rriT = \/p^+m% is the particle's transverse mass. 
The viscous deviation from local equilibrium is given by 

Hi 



Sfi(x,p) = /cq,j(p, x)(lTfcq,i(P' X )) 



P^p y TT^ix) 



2 T 2 e+p 



2T 2 (x) (e(x)+p{x)) 
(16) 



The approximation in the second line is not used in our 
numerical results but it holds (within the line thickness 
in all of our corresponding plots) since (lT/eq) deviates 
from 1 only when p<T where the last factor is small. 
With it the spectrum (fT2"l) takes the instructive form 



E- 



,d 3 N t 



d 3 p (27r) J 



p.d 3 a(x)f eq , i (x >P ){l+^^) 



+P< 
(17) 



The viscous correction is seen to be proportional to 
n 1 * 1 ' (x) on the freeze-out surface (normalized by the equi- 
librium enthalpy e+p) and to increase quadratically with 
the particle's momentum (normalized by the temperature 
T). At large pr, the viscous correction can exceed the 
equilibrium contribution, indicating a breakdown of vis- 
cous hydrodynamics. In that domain, particle spectra 
can not be reliably computed with viscous fluid dynam- 
ics. The limit of applicability depends on the actual value 
of tt^ /(e+p) and thus on the specific dynamical condi- 
tions encountered in the heavy-ion collision. 

The viscous correction to the spectrum in Eq. (JT7J) 
reads explicitly 

PuPu^" = mr(cosh 2 (?/-?7)7r rT + sinh 2 (y~77)r 2 7r ,?,? ) 
- 2uit cosh(y-rj) (p x ir TX + p y K Ty ) 

+ (pl^ x + 2p xPy ^y + P yy). (is) 

We can use the expressions given in Appendix 2 of 
Ref. [iH (in particular Eqs. (A22) in that paper) to re- 
express this in terms of the three independent compo- 
nents of TT mn for which we choose 



TT XX +ir VV , A = TT XX -TT yy . (19) 



This choice is motivated by our numerical finding (see 

! and ir vv are 



Fig. 2 in [3J and Sec. ITVC) that tt 7 ™, 
about an order of magnitude larger than all other com- 
ponents of 7r mn , and that in the azimuthally symmetric 
limit of central (b = 0) heavy-ion collisions the azimuthal 
average of A vanishes (see Eq. (|C2|>): (A)^ = 0. We find 
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P»P^ V = * m 



E 
A 



2/0 v.2/ \ i\ P T i / ,sin(0 p +^) p T sin(20 p ) 

777 T (2cosh (j/-r?)-l) - 2— m r cosh(y-7 ? ) . ^- + . . '\ 
x v± sm(2<p v ) uj_ sin(20 tI ) 

2 ,2/ \ „PT { sin((f)p+<f> v ) vj_sm(^ p -^) v ) 

771^ COSh U— 77 — 2 777 T COShlU— 77) . j- , . , . 

T Ky " v±_ V sin(2(7i t) ) 2 tan(20„) 



Ptttit cosh(y—r))vj_ 



sin 



bp-M Pt sin(2((/)p-^)) 



sin^c/)^) 



sin^c/)^) 
I 



2 «i V 2 / sm(2</>„) 



(20) 



Due to longitudinal boost-invariance, the integration 
over space-time rapidity 77 in Eq. ()12|) can be done an- 
alytically, resulting in a series of contributions involving 
modified Bessel functions [H, [Hj|. VISH2+1 docs not 
exploit this possibility and instead performs this and all 
other integrations for the spectra numerically. 

Once the spectrum (fT2"j) has been computed, a Fourier 
decomposition with respect to the azimuthal angle 4> p 
yields the anisotropic flow coefficients. For collisions be- 
tween equal spherical nuclei followed by longitudinally 
boost-invariant expansion of the collision fireball, only 
even-numbered coefficients contribute, the "elliptic flow" 
7J2 being the largest and most important one: 



d 3 Ni , s 



dN, 



d 3 p 
1 



dyprdpT 1 



■(b) 



(21) 



dN, 



2n dypxdpT 



[l + 2v 2 (p T ;b) cos(20 p ) 



In practice it is evaluated as the cos(20 p )-moment of the 
final particle spectrum, 



v 2 (j>t) = (cos(20 p )) 



/ d(j) p cos(20 ; 



dN 



P ' dypTdpTdtpp 



I' 



f, dN 

P dyp T dp T d(j> p 



, (22) 



where, according to Eq. (Ti"2")) . the particle spectrum is a 
sum of a local equilibrium and a non-equilibrium contri- 
bution (to be indicated symbolically as N = N cq + 5N). 



III. CENTRAL COLLISIONS 

A. Hydrodynamic evolution 

Even without transverse flow initially, the boost- 
invariant longitudinal expansion leads to a non- vanishing 
initial stress tensor a mn which generates non-zero target 
values for three components of the shear viscous pres- 



sure tensor: t 2 tt vv = 2 , n x 



J^. Inspection of 



the source terms in Eqs. then reveals that the ini- 

tially negative T 2 -K m reduces the longitudinal pressure, 
thus reducing the rate of cooling due to work done by 
the latter, while the initially positive values of n xx and 
w vy add to the transverse pressure and accelerate the de- 
volpment of transverse flow in x and y directions. As the 
fireball evolves, the stress tensor a mn receives additional 



contributions involving the transverse flow velocity and 
its derivatives (see Eq. (|A2j) ) which render an analytic 
discussion of its effects on the shear viscous pressure im- 
practical. 

Figure [5] shows what one gets numerically. Plotted are 
the source terms and (|5J), averaged over the trans- 
verse plane with the energy density as weight function, 
as a function of time, for evolution of central Cu+Cu 
collisions with two different equations of state, EOS I 
and SM-EOS Q. (In central collisions (|«S TX |) = (|«S T *|).) 
One sees that the initially strong viscous reduction of 
the (negative) source term S TT , which controls the cool- 
ing by longitudinal expansion, quickly disappears. This 
is due to a combination of effects: while the magni- 
tude of T 2 -K m decreases with time, its negative effects 
are further compensated by a growing positive contri- 
bution T^d x (pv x )+dy(pv y )) arising from the increasing 
transverse flow gradients. In contrast, the viscous in- 
crease of the (positive) transverse source term S TX per- 
sists much longer, until about 5fm/c. After that time, 
however, the viscous correction switches sign (clearly vis- 
ible in the upper inset in the right panel of Fig. [2t>) and 
turns negative, thus reducing the transverse acceleration 
at late times relative to the ideal fluid case. We can 
summarize these findings by stating that shear viscosity 
reduces longitudinal cooling mostly at early times while 
causing initially increased but later reduced acceleration 
in the transverse direction. Due to the general small- 
ness of the viscous pressure tensor components at late 
times, the last-mentioned effect (reduced acceleration) is 
not very strong. 

The phase transition in SM-EOS Q is seen to cause an 
interesting non-monotonic behaviour of the time evolu- 
tion of the source terms (right panel in Fig. [2]), leading 
to a transient increase of the viscous effects on the longi- 
tudinal source term while the system passes through the 
mixed phase. 

The viscous slowing of the cooling process at early 
times and the increased rate of cooling at later times 
due to accelerated transverse expansion are shown in Fig- 
urc[3] The upper set of curves shows what happens in the 
center of the fireball. For comparison we also show curves 
for boost-invariant longitudinal Bjorkcn expansion with- 
out transverse flow, labeled "(0+l)-d hydro". These are 
obtained with flat initial density profiles for the same 
value eo (no transverse gradients). The dotted green 




FIG. 2: (Color online) Time evolution of the hydrodynamic source terms (|7I9[) . averaged over the transverse plane, for central 
Cu+Cu collisions, calculated with EOS I in the left panel and with SM-EOS Q in the right panel. The smaller insets blow 
up the vertical scale to show more detail. The dashed blue lines are for ideal hydrodynamics with eo = 30 GeV/fm 3 and 
To = 0.6fm/c. Solid red lines show results from viscous hydrodynamics with identical initial conditions and 2 = j- ~ 0.08, 
tw = Jj? ss 0.24 ( 200 Mcv ) fm/c. The positive source terms drive the transverse expansion while the negative ones affect the 
longitudinal expansion. 




FIG. 3: (Color online) Time evolution of the local temperature in central Cu+Cu collisions, calculated with EOS I (left) and 
SM-EOS Q (right), for the center of the fireball (r — 0, upper set of curves) and a point near the edge (r = 9fm, lower set of 
curves). Same parameters as in Fig. [5] See text for discussion. 



line in the left panel shows the well-known T^r -1 / 3 
behaviour of the Bjorken solution of ideal fluid dynam- 
ics [HI , modified in the right panel by the quark-hadron 
phase transition where the temperature stays constant in 
the mixed phase. The dash-dotted purple line shows the 
slower cooling in the viscous (0+l)-dimensional case [l7l |. 
due to reduced work done by the longitudinal pressure. 
The expansion is still boost-invariant a la Bjorken [54[ 
(as it is for all other cases discussed in this paper), but 
viscous effects generate entropy, thereby keeping the tem- 
perature at all times higher than for the adiabatic case. 
The dashed blue (ideal) and solid red (viscous) lines for 
the azimuthally symmetric (l+l)-dimcnsional case show 



the additional cooling caused by transverse expansion. 
Again the cooling is initially slower in the viscous case 
(solid red), but at later times, due to faster build-up of 
transverse flow by the viscously increased transverse pres- 
sure, the viscous expansion is seen to cool the fireball 
center faster than ideal hydrodynamics. (Note also the 
drastic reduction of the lifetime of the mixed phase by 
transverse expansion; due to increased transverse flow 
and continued acceleration in the mixed phase from vis- 
cous pressure gradients, it is even more dramatic in the 
viscous than the ideal case.) The curves for r = 9 fm cor- 
roborate this, showing that the temperature initially in- 
creases with time due to hot matter being pushed from 




FIG. 4: (Color online) Time evolution of the local entropy density for central Cu+Cu collisions, calculated with EOS I (left) 
and SM-EOS Q (right), for the center of the fireball (r = 0, upper set of curves) and a point at r = 3fm (lower set of curves). 
Same parameters and color coding as in Fig. [3] See text for discussion. 



the center towards the edge, and that this temperature 
increase happens more rapidly in the viscous fluid (solid 
red lines), due to the faster outward transport of matter 
in this case. 

Figure [4] shows how the features seen in Fig. [3] manifest 
themselves in the evolution of the entropy density. (In 
the QGP phase s^T 3 .) The double-logarithmic presen- 
tation emphasizes the effects of viscosity and transverse 
expansion on the power law s(t) ~ r~ a : One sees that 
the r _1 scaling of the ideal Bjorken solution is flattened 
by viscous effects, but steepened by transverse expan- 
sion. As is well-know, it takes a while (here about 3fm/c) 
until the transverse rarefaction wave reaches the fireball 
center and turns the initially 1-dimcnsional longitudinal 
expansion into a genuinely 3-dimensional one. When this 
happens, the power law s(r) ~ r~ a changes from a = 1 in 
the ideal fluid case to a > 3 [l| . Here 3 is the dimension- 
ality of space, and the fact that a becomes larger than 
3 reflects rclativistic Lorentz-contraction effects through 
the transverse-flow-rclated 7j_-factor that keeps increas- 
ing even at late times. In the viscous case, a changes 
from 1 to 3 sooner than for the ideal fluid, due to the 
faster growth of transverse flow. At late times the s(r) 
curves for ideal and viscous hydrodynamics are almost 
perfectly parallel, indicating that very little entropy is 
produced during this late stage. 

In Figure [5] we plot the evolution of temperature in 
r— t space, in the form of constant-T surfaces. Again 
the two panels compare the evolution with EOS I (left) 
to the one with SM-EOS Q (right). In the two halves 
of each panel we directly contrast viscous and ideal fluid 
evolution. (The light gray lines in the right halves are re- 
flections of the viscous temperature contours in the left 
halves, to facilitate comparison of viscous and ideal fluid 
dynamics.) Beyond the already noted fact that at r = Q 
the viscous fluid cools initially more slowly (thereby giv- 
ing somewhat longer life to the QGP phase) but later 



more rapidly (thereby freezing out earlier), this figure 
also exhibits two other noteworthy features: (i) Moving 
from r = outward, one notes that contours of larger ra- 
dial flow velocity are reached sooner in the viscous than 
in the ideal fluid case; this shows that radial flow builds 
up more quickly in the viscous fluid. This is illustrated 
more explicitly in Fig. [5] which shows the time evolution 
of the radial velocity (v±), calculated as an average over 
the transverse plane with the Lorentz contracted energy 
density 7j_e as weight function, (ii) Comparing the two 
sets of temperature contours shown in the right panel of 
Fig. [51 one sees that viscous effects tend to smoothen any 
structures related to the (first order) phase transition in 
SM-EOS Q. The reason for this is that, with the dis- 
continuous change of the speed of sound at either end of 
the mixed phase, the radial flow velocity profile develops 
dramatic structures at the QGP-MP and MP-HRG inter- 
faces [ii]]. This leads to large velocity gradients across 
these interfaces (as can be seen in the right panel of Fig. [5] 
in its lower right corner which shows rather twisted con- 
tours of constant radial flow velocity) , inducing large vis- 
cous pressures which drive to reduce these gradients (as 
seen in lower left corner of that panel). In effect, shear 
viscosity softens the first-order phase transition into a 
smooth but rapid cross-over transition. 



These same viscous pressure gradients cause the fluid 
to accelerate even in the mixed phase where all thermo- 
dynamic pressure gradients vanish (and where the ideal 
fluid therefore docs not generate additional flow). As a 
result, the lifetime of the mixed phase is shorter in vis- 
cous hydrodynamics, as also seen in the right panel of 
Figure 



10 



Cu+Cu, b=0, EOS I Cu+Cu, b=0, SM-EOS Q 




-5 5 -5 5 

r(fm) r ( fm ) 



FIG. 5: (Color online) Surfaces of constant temperature T and constant transverse flow velocity v± for central Cu+Cu collisions, 
evolved with EOS I (left panel) and SM-EOS Q (right panel). In each panel, results from viscous hydrodynamics in the left 
half are directly compared with the corresponding ideal fluid evolution in the right half. (The thin isotherm contours in the 
right halves of each panel are reflected from the left halves, for easier comparison.) The lines of constant v± are spaced by 
intervals of 0.1, from the inside outward, as indicated by the numbers near the top of the figures. The right panel contains 
two isotherms for T c = 164 MeV, one separating the mixed phase (MP) from the QGP at energy density eQ = 1.6 GeV/fm 3 , the 
other separating it from the hadron resonance gas (HRG) at energy density eH = 0.45 GeV/fm 3 . See text for discussion. 
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FIG. 6: (Color online) Time evolution of the average radial flow velocity (vt) = (v±) in central Cu+Cu collisions, calculated 
with EOS I (left panel) and SM-EOS Q (right panel). Solid (dashed) lines show results from ideal (viscous) fluid dynamics. 
The initially faster rate of increase reflects large positive shear viscous pressure in the transverse direction at early times. The 
similar rates of increase at late times indicate the gradual disappearance of shear viscous effects. In the right panel, the curves 
exhibit a plateau from 2 to 4fm/c, reflecting the softening of the EOS in the mixed phase. 



B. Final particle spectra 



After obtaining the freeze-out surface, we calculate the 
particle spectra from the generalized Cooper-Frye for- 
mula (12]), using the AZHYDRO algorithm [U for the 
integration over the freeze-out surface S. For calcula- 
tions with EOS I which lacks the transition from mass- 
less partons to hadrons, we cannot compute any hadron 
spectra. For illustration we instead compute the spectra 



of hypothetical massless bosons ("gluons"). They can be 
compared with the pion spectra from SM-EOS Q which 
can also, to good approximation, be considered as mass- 
less bosons. 

The larger radial flow generated in viscous hydrody- 
namics, for a fixed set of initial conditions, leads, of 
course, to flatter transverse momentum spectra [30l HH, 
HH (at least at low px where the viscous correction 5fi 
to the distribution function can be neglected in (TT!?|) ). 
This is seen in Figure [JJ by comparing the dotted and 
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FIG. 7: (Color online) Mid-rapidity particle spectra for central Cu+Cu collisions, calculated with EOS I (left, gluons) and 
with SM-EOS Q (right, 7r~, K + and p). The solid blue (red dashed) lines are from ideal (viscous) hydrodynamics. The purple 
dotted lines show viscous hydrodynamic spectra that neglect the viscous correction Sfi to the distribution function in Eq. (|12|l , 
i.e. include only the effects from the larger radial flow generated in viscous hydrodynamics. 



solid lines. This comparison also shows that the viscous 
spectra lie systematically above the ideal ones, indicating 
larger final total multiplicity. This reflects the creation 
of entropy during the viscous hydrodynamic evolution. 
As pointed out in [33l . |34| . this requires a retuning of 
initial conditions (starting the hydrodynamic evolution 
later with smaller initial energy density) if one desires to 
fit a given set of experimental p^-spectra. Since we here 
concentrate on investigating the origins and detailed me- 
chanics of viscous effects in relativistic hydrodynamics, 
we will not explore any variations of initial conditions. 
All comparisons between ideal and viscous hydrodynam- 
ics presented here will use identical starting times To and 
initial peak energy densities eo- 

The viscous correction Sfi in Eqs. (|12I16[) depends on 
the signs and magnitudes of the various viscous pressure 
tensor components along the freeze-out surface, weighted 
by the equilibrium part / eqi i of the distribution function. 
Its effect on the final py-spectra (even its sign!) is not a 
priori obvious. Teaney [lq |. using a blast- wave model 
to evaluate the velocity stress tensor = n^ v /(2r)), 
found that the correction is positive, growing quadrat- 
ically with px- Romatschke et al. [3J, [36[ did not 
break out separately the contributions from larger ra- 
dial flow in / cq i and from Sfi. Dusling and Teaney [38| . 
solving a slightly different set of viscous hydrodynamic 
equations and using a different (kinetic) freeze-out cri- 
terium to determine their decoupling surface, found a 
(small) positive effect from Sfi on the final pion spec- 
tra, at least up to pr = 2GeV /c, for freeze-out around 
Tdcc ~ 130 MeV, turning weakly negative when their ef- 
fective freeze-out temperature was lowered to below 100 
MeV. The dashed lines in Figure [7] show that in our cal- 
culations for pT^2GeV/c the effects from Sfi have an 
overall negative sign, leading to a reduction of the px- 
spectra at large pr relative to both the viscous spectra 
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FIG. 8: (Color online) Ratio of the viscous correction SN, 
resulting from the non-equilibrium correction 5f, Eq. (|16[1 . 
to the distribution function at freeze-out, to the equilibrium 
spectrum N e q = dN eq /(dyd 2 pT) calculated from Eq. (I12|l by 
setting Sf = 0. The gluon curves are for evolution with EOS I, 
the curves for n~ , K + and p are from calculations with SM- 
EOS Q. 

without Sfi and the ideal hydrodynamic spectra. This is 
true for all particle species, irrespective of the EOS used 
to evolve the fluid. 

It turns out that, when evaluating the viscous cor- 
rection Sf in Eq. (|16|) with the help of Eq. (|2T)]) . large 
cancellations occur between the first and second line in 
Eq. (f20|) . [After azimuthal integration, the contribution 
to Sf from the third line ~ A vanishes identically for cen- 
tral collisions.] These cancellations cause the final result 
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to be quite sensitive to small numerical errors in the cal- 
culation of T 2 7r w and £ = tt xx +it vv . Increased numerical 



stability is achieved by trading t 2 tt vv for ir TT = T 2 7r''''+S 
and using instead of Eq. ([20]) the following expression: 



+ A 



rn3p{2 cosh 2 (y— rj)— l) — 2— ror cosh(y— rf) 
— m 2 " sinh (y— rj) + prmx coah(y—rj)vj_ 
Ptttit cosh(y— rj)v 



sm((j) p +(j) v ) p\ sin(20 p ) 



sin(2^ v ) Vj_ sin(2</> t ,) 
sa\((j) p -cf> v ) jfy / _ sin(20 p ) 
tan(20„) 2 \ sin(2^>„) 



sin(^p-0„) p| sin(2(0 p -0„)) 



sin(20„) 



sin(2</>„) 



r 



(23) 



The first and second lines of this expression are now much 
smaller than before and closer in magnitude to the final 
net result for p tl p v TT fJ/l/ . This improvement carries over to 
non-central collisions as discussed in Sec. IIVD1 where we 
also show the individual contributions from ir TT , E and 
A to the particle spectra. To be able to use Eq. ([23]) . 
the numerical code should directly evolve, in addition 
to 7r TT , tt tx , and 7r ry which are needed for the velocity 
finding algorithm (see Appendix IB")) . the components ir xx 
and it vv . Otherwise these last two components must be 
computed from the evolved ir mn components using the 
transversality and tracelessness constraints which neces- 
sarily involves the amplification of numerical errors by 
division with small velocity components. 

In Figure [8] we explore the non-equilibrium contribu- 
tion to the final hadron spectra in greater detail. The 
figure shows that the non-equilibrium effects from 5fi 
are largest for massless particles and, at high px, de- 
crease in magnitude with increasing particle mass. The 
assumption |<5/|<C/eqj which underlies the viscous hy- 
drodynamic formalism, is seen to break down at high px, 
but to do so later for heavier hadrons than for lighter 
ones. Once the correction exceeds 0(50%) (indicated by 
the horizontal dashed line in Fig. [5]), the calculated spec- 
tra can no longer be trusted. 

In contrast to viscous hydrodynamics, ideal fluid dy- 
namics has no intrinsic characteristic that will tell us 
when it starts to break down. Comparison of the calcu- 
lated elliptic flow «2 from ideal fluid dynamics with the 
experimental data from RIIIC [l[ suggests that the ideal 
fluid picture begins to break down above pr — 1.5 GeV/c 
for pions and above pt — 2 GeV/c for protons. This phe- 
nomenological hierarchy of thresholds where viscous ef- 
fects appear to become essential is qualitatively consis- 
tent with the mass hierarchy from viscous hydrodynamics 
shown in Fig. [8l 

In the region <pr ^$ 1-5 GeV/c, the interplay be- 
tween rriT- and p^-dependent terms in Eq. (|20[) is subtle, 
causing sign changes of the viscous spectral correction 
depending on hadron mass and pr (see inset in Fig. [5J . 
The fragility of the sign of the effect is also obvious from 
Fig. 8 in the work by Dusling and Teaney [38| where it 



is shown that in this px region the viscous correction 
changes sign from positive to negative when freeze-out is 
shifted from earlier to later times (higher to lower freeze- 
out temperature). Overall, we agree with them that the 
viscous correction effects on the py-spectra are weak in 
this region [3a |. We will see below that a similar state- 
ment does not hold for the elliptic flow. 



IV. NON-CENTRAL COLLISIONS 
A. Hydrodynamic evolution 

We now take full advantage of the ability of VISH2+1 
to solve the transverse expansion in 2 spatial dimensions 
to explore the anisotropic fireball evolution in non-central 
heavy-ion collisions. Similar to Fig. [S] for 6 = 0, Figure [5] 
shows surfaces of constant temperature and radial flow 
for Cu+Cu collisions at 6= 7fm, for the equation of state 
SM-EOS Q. The plots show the different evolution into 
and perpendicular to the reaction plane and compare 
ideal with viscous fluid dynamics. Again, even a minimal 
amount of shear viscosity ( f = ^ ) is seen to dramatically 
smoothen all structures related to the existence of a first- 
order phase transition in the EOS. However, in distinc- 
tion to the case of central collisions, radial flow builds 
up at a weaker rate in the peripheral collisions and never 
becomes strong enough to cause faster central cooling at 
late times than seen in ideal fluid dynamics (bottom row 
in Fig. [9]). The viscous fireball cools more slowly than 
the ideal one at all times and positions, lengthening in 
particular the lifetime of the QGP phase, and it grows 
to larger transverse size at freeze-out. [Note that this 
does not imply larger transverse HBT radii than for ideal 
hydrodynamics (something that -in view of the "RHIC 
HBT Puzzle" [l|- would be highly desirable): the larger 
geometric size is counteracted by larger radial flow such 
that the geometric growth, in fact, does not lead to larger 
transverse homogeneity lengths [34j.] 

While Figure [5] gives an impression of the anisotropy 
of the fireball in coordinate space, we study now in 
Fig. [TO] the evolution of the flow anisotropy (|Uc| — \ v y\)- 
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FIG. 9: (Color online) Surfaces of constant temperature T and constant transverse flow velocity v± for semi-peripheral Cu+Cu 
collisions at 6 = 7fm, evolved with SM-EOS Q. In the top row we contrast ideal (left panel) and viscous (right panel) fluid 
dynamics, with a cut along the x axis (in the reaction plane) shown in the right half while the left half shows a cut along the y 
axis (perpendicular to the reaction plane). In the bottom row we compare ideal and viscous evolution in the same panel, with 
cuts along the x (y) direction shown in the left (right) panel. See Fig. [5] for comparison with central Cu+Cu collisions. 



In central collisions this quantity vanishes. In ideal hy- 
drodynamics it is driven by the anisotropic gradients of 
the thermodynamic pressure. In viscous fluid dynamics, 
the source terms (|8I9[) . whose difference is shown in the 
bottom row of Fig. [TUJ receive additional contributions 
from gradients of the viscous pressure tensor which con- 
tribute their own anisotropies. Fig. llOl demonstrates that 
these additional anisotropies increase the driving force for 
anisotropic flow at very early times (t— tq < 1 fm/c), but 
reduce this driving force throughout the later evolution. 
At times r— Tn>2fm/c the anisotropy of the effective 
transverse pressure even changes sign and turns nega- 
tive, working to decrease the flow anisotropy. As a con- 
sequence of this, the buildup of the flow anisotropy stalls 
at r— To ~ 2.5 fm/c (even earlier for SM-EOS Q where 
the flow buildup stops as soon as the fireball medium en- 
ters the mixed phase) and proceeds to slightly decrease 



thcraftcr. This happens during the crucial period where 
ideal fluid dynamics still shows strong growth of the flow 
anisotropy. By the time the fireball matter decouples, 
the average flow velocity anisotropy of viscous hydro lags 
about 20-25% behind the value reached during ideal fluid 
dynamical evolution. 

These features are mirrored in the time evolution of 
the spatial eccentricity e x = r^q—sj (calculated by aver- 
aging over the transverse plane with the energy density 
e(x) as weight function [4J] and shown in the top row 
of Fig. Hip and of the momentum anisotropies e p and e' p 
(shown in the bottom row). The momentum anisotropy 

/rpXX_rpyy\ I ■ 

(p = T °„i T in [55fl measures the anisotropy of the trans- 
verse momentum density due to anisotropies in the col- 
lective flow pattern, as shown in top row of Fig. 1101 it in- 
cludes only the ideal fluid part of the energy momentum 
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FIG. 10: (Color online) Time evolution of the transverse flow anisotropy (|i>a:| — |i>ji|) (top row) and of the anisotropy in the 
transverse source term (|<S Ta: | — |<S Ta |) (bottom row). Both quantities are averaged over the transverse plane, with the Lorentz- 
contracted energy density j± as weight function. The left (right) column shows results for EOS I (SM-EOS Q), with solid 
(dashed) lines representing ideal (viscous) fluid dynamical evolution. 



tensor. The total momentum anisotropy e' p = | Txx+Tyi) | , 
similarly defined in terms of the total energy momentum 
tensor ~Tq V +tt^ v , additionally counts anisotropic 
momentum contributions arising from the viscous pres- 
sure tensor. Since the latter quantity includes effects 
due to the deviation Sf of the local distribution func- 
tion from its thermal equilibrium form which, accord- 
ing to Eq. (1121) . also affects the final hadron momen- 
tum spectrum and elliptic flow, it is this total momen- 
tum anisotropy that should be studied in viscous hydro- 
dynamics if one wants to understand the evolution of 
hadron elliptic flow. In other words, in viscous hydrody- 
namics hadron elliptic flow is not simply a measure for 
anisotropies in the collective flow velocity pattern, but 
additionally reflects anisotropies in the local rest frame 
momentum distributions, arising from deviations of the 
local momentum distribution from thermal equilibrium 
and thus being related to the viscous pressure. 

Figure [IT] correlates the decrease in time of the spa- 
tial eccentricity e x with the buildup of the momentum 
anisotropies e p and e' . In viscous dynamics the spatial 
eccentricity is seen to decrease initially faster than for 
ideal fluids. This is less a consequence of anisotropies in 
the large viscous transverse pressure gradients at early 
times than due to the faster radial expansion caused by 
their large overall magnitude. In fact, it was found a 



while ago [56[ that for a system of free-streaming partons 
the spatial eccentricity falls even faster than the viscous 
hydrodynamic curves (solid lines) in the upper row of Fig- 
ure[llj The effects of early pressure gradient anisotropics 
is reflected in the initial growth rate of the flow-induced 
momentum anisotropy e p which is seen to slightly ex- 
ceed that observed in the ideal fluid at times up to about 

1 fm/c after the beginning of the transverse expansion 
(bottom panels in Fig. [TTj) . This parallels the slightly 
faster initial rise of the flow velocity anisotropy seen in 
the top panels of Fig. [TUJ Figure [TU] also shows that in 
the viscous fluid the flow velocity anisotropy stalls about 

2 fm/c after start and remains about 25% below the fi- 
nal value reached in ideal fluid dynamics. This causes 
the spatial eccentricity of the viscous fireball to decrease 
more slowly at later times than that of the ideal fluid 
(top panels in Fig. [TT|) which, at late times, features a 
significantly larger difference between the horizontal (x) 
and vertical (y) expansion velocities. 

It is very instructive to compare the behaviour of 
the flow-induced ideal-fluid contribution to the momen- 
tum anisotropy, e p , with that of the total momentum 
anisotropy e' p . At early times they are very different, 
with e' p being much smaller than e p and even turning 
slightly negative at very early times (see insets in the 
lower panels of Fig. [11]). This reflects very large nega- 
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FIG. 11: (Color online) Time evolution for the spatial eccentricity e x , momentum anisotropy e p and total momentum anisotropy 
e' p (see text for definitions), calculated for fe = 7fm Cu+Cu collisions with EOS I (left column) and SM-EOS Q (right column). 
Dashed lines are for ideal hydrodynamics while the solid and dotted lines show results from viscous hydrodynamics. See text 
for discussion. 



tive contributions to the anisotropy of the total energy 
momentum tensor from the shear viscous pressure whose 
gradients along the out-of-plane direction y strongly ex- 
ceed those within the reaction plane along the x direc- 
tion. At early times this effect almost compensates for 
the larger in-plane gradient of the thermal pressure. The 
negative viscous pressure gradient anisotropy is responsi- 
ble for reducing the growth of flow anisotropics, thereby 
causing the flow-induced momentum anisotropy e p to sig- 
nificantly lag behind its ideal fluid value at later times. 
The negative viscous pressure anisotropics responsible for 
the difference between e p and e' p slowly disappear at later 
times, since all viscous pressure components then become 
very small (see Fig. [131 below). 

The net result of this interplay is a total momentum 
anisotropy in Cu+Cu collisions (i.e. a source of elliptic 
flow V2) that for a "minimally" viscous fluid with 2 = J- 
is 40-50% lower than for an ideal fluid, at all except the 
earliest times (where it is even smaller) . The origin of this 
reduction changes with time: Initially it is dominated by 
strong momentum anisotropics in the local rest frame, 
with momenta pointing preferentially out-of-plane, in- 
duced by deviations from local thermal equilibrium and 
associated with large shear viscous pressure. At later 
times, the action of these anisotropic viscous pressure 



gradients integrates to an overall reduction in collective 
flow anisotropy, while the viscous pressure itself becomes 
small; at this stage, the reduction of the total momentum 
anisotropy is indeed mostly due to a reduced anisotropy 
in the collective flow pattern while momentum isotropy 
in the local fluid rest frame is approximately restored. 

B. Elliptic flow V2 of final particle spectra 

The effect of the viscous suppression of the total mo- 
mentum anisotropy e' on the final particle elliptic flow 
is shown in Figure [T2l Even for the "minimal" viscosity 
2 = t~ considered here one sees a very strong suppression 
of the differential elliptic flow V2{pt) from viscous evo- 
lution (dashed lines) compared to the ideal fluid (solid 
lines). Both the viscous reduction of the collective flow 
anisotropy (whose effect on vi is shown as the dotted 
lines) and the viscous contributions to the anisotropy of 
the local momentum distribution (embodied in the term 
Sf in Eq. (fT2"|) ) play big parts in this reduction. The runs 
with EOS 1 (which is a very hard EOS) decouple more 
quickly than those with SM-EOS Q; correspondingly, the 
viscous pressure components are still larger at freeze-out 
and the viscous corrections Sf to the distribution func- 
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FIG. 12: (Color online) Differential elliptic flow V2(pr) for Cu+Cu collisions at 6 = 7fm. Left panel: Gluons from evolution 
with EOS I. Right panel: ty~ , K + , and p from evolution with SM-EOS Q. Dashed lines: ideal hydrodynamics. Solid lines: 
viscous hydrodynamics. Dotted lines: viscous hydrodynamics without non-equilibrium distortion 8f of distribution function at 
freeze- out. 
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FIG. 13: (Color online) Left panel: Time evolution of the various components of the shear viscous pressure tensor, normalized 
by the enthalpy and averaged in the transverse plane over the thermalized region inside the freeze-out surface [5{j. Note that 
the normalizing factor e+p ~ T 4 decreases rapidly with time. Right panel: Comparison of the full viscous hydrodynamic source 
terms S mn , averaged over the transverse plane, with their approximations given in Eqs. (OH, as a function of time. The thin 
gray lines indicate the corresponding source terms in ideal fluid dynamics. 



tion play a bigger role. With SM-EOS Q the fireball 
doesn't freeze out until ir mn has become very small (see 
Fig.[T3"lbelow), resulting in much smaller corrections from 
Sf (difference between dashed and dotted lines in Fig.[T2"|) 
[57l | . On the other hand, due to the longer fireball life- 
time the negatively anisotropic viscous pressure has more 
time to decelerate the buildup of anisotropic flow, so vi 
is strongly reduced because of the much smaller flow- 
induced momentum anisotropy e p . 

The net effect of all this is that, for Cu+Cu collisions 
and in the soft momentum region pt < 1.5 GeV/c, the 
viscous evolution with 2 = -L leads to a suppression of 
i>2 by about a factor 2 (58j . in both the slope of its px- 
dependence and its pr-integrated value. (Due to the flat- 



ter pT-spectra from the viscous dynamics, the effect in 
the pT-integrated v-i is not quite as large as for V2{pt) at 
fixed Pt-) 



C. Time evolution of the viscous pressure tensor 
components and hydrodynamic source terms 

In Figure [13] we analyze the time evolution of the vis- 
cous pressure tensor components and the viscous hydro- 
dynamic source terms on the r.h.s. of Eqs. ([6|). As 
already mentioned, the largest components of 7r m " are 
t 2 tt vv , ir xx and ir yy (see Fig. 2 in [33| and left panel of 
Fig. [TJ [H). At early times, both r 2 ^ and the sum 
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E = n xx +ir vv reach (with opposite signs) almost 20% of 
the equilibrium enthalpy e+p. At this stage all other 
components of n are at least an order of magnitude 
smaller (see inset). The largest of these small compo- 
nents is the difference A = ir xx — ir vv which we choose 
as the variable describing the anisotropy of the viscous 
pressure tensor in non-central collisions. At late times 
(t— To > 5 fm/c), when the large components of 7r m ™ have 
strongly decreased, A becomes comparable to them in 
magnitude. As a fraction of the thermal equilibrium en- 
thalpy e+p^T 4 which sets the scale in ideal fluid dy- 
namics and which itself decreases rapidly with time, all 
viscous pressure components are seen to decrease with 
time. In a fluid with a set ratio rj/s, viscous effects thus 
become less important with time. In real life, however, 
the ratio rj/s depends itself on temperature and rises dra- 
matically during the quark-hadron phase transition and 
below 0, EH. Shear viscous effects will therefore be 
larger at late times than considered here. The conse- 
quences of this will be explored elsewhere. 

The observation that many components of 7r Tnn are 
very small throughout the fireball evolution underlies the 
validity of the approximation of the hydrodynamic source 
terms given in the second lines of Eqs. (JTHH]). The excel- 
lent quality of this approximation is illustrated in the 
right panel of Fig. IT51 



D. Viscous corrections to final pion spectra and 
elliptic flow 

The large viscous reduction of the elliptic flow seen in 
Fig. [T2] warrants a more detailed analysis of the viscous 
corrections to the particle spectra and V2- In Fig. [TJ] 
we show, for Cu+Cu collisions at b = 7 fm evolved with 
SM-EOS Q, the time evolution of the independent com- 
ponents 7r TT , E and A of the viscous pressure tensor 7r m ™, 
normalized by the equilibrium enthalpy e+p, along the 
Idee =130 Me V decoupling surface plotted in the upper 
right panel of Fig. [9l Solid (dashed) lines show the be- 
haviour along the x (y) direction (right (left) half of the 
upper right panel in Fig. [5]). We see that generically all 
three of these viscous pressure components are of similar 
magnitude, except for E which strongly dominates over 
the other two during the first 2 fm/c after the beginning 
of the expansion stage. However, since most particle pro- 
duction, especially that of 1ow-pt particles, occurs at late 
times (r>4fm/c for 6 = 7fm/c Cu+Cu, see Fig. [5] and 
the discussion around Fig. 27 in Ref. [l(), the regions 
where E is large do not contribute much. As far as the 
non-equilibrium contribution to the spectra is concerned, 
we can thus say that the viscous pressure at freeze-out 
is of the order of a few percent of e+p. The anisotropy 
term A is even smaller, due to cancellations between the 
in-plane (x) and out-of-plane (y) contributions when in- 
tegrating over the azimuthal angle in Eq. (TT2")) . 

These viscous pressure components generate the non- 
equilibrium contribution 5f to the distribution func- 




FIG. 14: Time evolution of n TT , E, and A, as well as that 
of the transverse velocity, along the decoupling surface for 
b= 7fm Cu+Cu collisions in viscous hydrodynamics with SM- 
EOS Q, as shown in Fig. Solid (dashed) lines represent cuts 
along the x (</> = 0) and y (<^> = §) directions. 



tion on the freeze-out surface according to Eqs. (Til)]) 
and (|23p . resulting in a corresponding viscous correc- 
tion to the azimuthally integrated particle spectrum 

SN= fdcf) p s(E^y Figure [15] shows these non- 
equilibrium contributions for pions, normalized by the az- 
imuthally averaged equilibrium part N cq = J d(f> p rf J^ q . 
We show both the total viscous correction and the indi- 
vidual contributions arising from the three independent 
pressure tensor components used in Eq. (|23|) and shown 
in Fig. M 

In the viscous correction, the term (|23[) (normalized 
by T 2 {e+p)) is weighted by particle production via the 
equilibrium distribution function f e<l (x,p). It is well- 
known (see Fig. 27 in Ref. [IJ) that for low-pr particles 
this weight is concentrated along the relatively flat top 
part of the decoupling surface in Fig. [51 corresponding to 
t> 5— 6 fm/c in Fig. 1141 In this momentum range, the 
contributions from ir TT , E and A to 5N/N eq are of simi- 
lar magnitude and alternating signs (see Fig. ll5[) . making 
the sign of the overall viscous correction to the spectra 
hard to predict. 

High-pT particles, on the other hand, come from those 
regions in the fireball which feature the largest transverse 
flow velocity at freeze-out. Fig.[5]shows that this restricts 
their emission mostly to the time interval 3 < t < 6 fm/c. 
In this region tt tt is negative, see Fig. [TU A detailed 
study of the different terms in Eq. (|23[) reveals that (after 
azimuthal integration) the expression multiplying ir TT is 
positive, hence the negative sign of 7r TT explains its nega- 
tive contribution to 5N/N eq at high pt, as seen in Fig. [TBI 
Figure [TBI also shows that in the region pt ^ 1 GcV/c the 
first line ~ ir TT in Eq. (j2li|) completely dominates the vis- 
cous correction to the spectra. We found that this in- 
volves additional cancellations between terms of oppo- 
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FIG. 15: (Color online) Viscous corrections to the azimuthally 
averaged pion spectrum resulting from individual components 
of the viscous pressure tensor n mn as indicated, as well as the 
total correction SN/N eq , for Cu+Cu collisions at b= 7fm with 
SM-EOS Q. The horizontal dashed line at -50% indicates the 
limit of validity. 



site sign (after azimuthal integration) inside the square 
brackets multiplying £ and A in the second and third 
line of Eq. (T25|) . Furthermore, the term ~ ir TT is the only 
contribution whose magnitude grows quadratically with 
Pt- For the contributions involving E and A, the appar- 
ent quadratic momentum dependence seen in Eq. (1231) is 
tempered by the integrations over space-time rapidity i] 
and azimuthal angle <f> in (|12p . resulting in only linear 
growth at large pr- 

In the absence of higher-order momentum anisotropies 
v n , n > 2, the elliptic flow V2{pr) can be easily computed 
from the momentum spectra in x (</> p =0) and y (</> P =f ) 
directions: 



2v 2 (pi 



N X — Ny 

N 

(iV^eq-A^eq) + (SN x -SN y 

N eq + 6N 



(24) 



where N = N eq +SN is shorthand for the azimuthally av- 
eraged spectrum dN / (2ir dy pxdpx), and N x y denote the 
Pt spectra along the x and y directions, respectively: 
N x = N XiCq +SN x = dyp f d p T dtt>p (<f> P =0), and similarly for 
N y with p =f . Equation (|24|) shows that V2 receives con- 
tributions from anisotropies in the equilibrium part of the 
distribution function / eq , which reflect the hydrodynamic 
flow anisotropy along the freeze-out surface, and from the 
viscous correction 8f, which reflects non-equilibrium mo- 
mentum anisotropies in the local fluid rest frame. The 
dashed line in Figure [16] shows the relative magnitude 
of these two anisotropy contributions, A f jVj: _^ V " , and 
compares it with the relative magnitude 4P- of the non- 
equilibrium and equilibrium contributions to the total, 
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FIG. 16: (Color online) Ratio of non-equilibrium and equilib- 
rium contributions to particle production (solid line) and to 
its momentum anisotropy (dashed line) , as a function of pr 
for pions from Cu+Cu collisions at 6 = 7fm with SM-EOS Q. 



(^-integrated pion spectrum for Cu+Cu at 6 = 7fm. We 
see that the non-equilibrium contribution to the momen- 
tum anisotropy V2 is always negative and larger in rel- 
ative magnitude than the non-equilibrium contribution 
to the azimuthally averaged spectrum. Since i>2 is a 
small quantity reflecting the anisotropic distortion of the 
single-particle spectrum, it reacts more sensitively than 
the spectrum itself to the (anisotropic) non-equilibrium 
contributions caused by the small viscous pressure ir mn 
on the decoupling surface. Furthermore, the viscous cor- 
rections to the 0p-integrated spectrum change sign as a 
function of pt, the corrections to i>2 are negative every- 
where, decreasing V2{pt) at all values of pr, but espe- 
cially at large transverse momenta. 



V. SENSITIVITY TO INPUT PARAMETERS 
AND LIMITS OF APPLICABILITY 

A. Initialization of n mn 

Lacking input from a microscopic model of the pre- 
cquilibrium stage preceding the (viscous) hydrodynamic 
one, one must supply initial conditions for the energy 
momentum tensor, including the viscous pressure ir mn . 
The most popular choice has been to initialize ir mn with 
its Navier-Stokes value, i.e. to set initially 7r m " = 2rja mn . 
Up to this point, this has also been our choice in the 
present paper. Ref. [36| advocated the choice Ti mn = 
at time To in order to minimize viscous effects and thus 
obtain an upper limit on rj/s by comparison with exper- 
imental data. In the present subsection we explore the 
sensitivity of the final spectra and elliptic flow to these 
different choices of initialization, keeping all other model 
parameters unchanged. 

Figure [T71 shows the time evolution of the viscous pres- 
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FIG. 17: (Color online) Similar to Figure 1131 but now comparing runs with different initial conditions. The thick lines 
reproduce the results from Figure [131 obtained with n mn = 2rja mn a,t initial time To, while thin lines of the same type show 
the corresponding results obtained by setting initially 7r mn = 0. The right panel shows the full viscous source terms, without 
approximation: {\S TX \) (dashed), {\S Ty \) (dotted), and (5 rT ) (dash-dotted). 



sure tensor and viscous hydrodynamic source terms for 
the two different initializations. Differences with respect 
to the results shown Fig. [T3] (which are reproduced in 
Fig. [IT] for comparison) are visible only at early times 
t — tq < 5tv pa 1 fm/c. After tv ~ 0.2 fm/c, the initial dif- 
ference ir mn —2r](j mn has decreased by roughly a factor 
1/e, and after several kinetic scattering times tv the 
hydrodynamic evolution has apparently lost all memory 
how the viscous terms were initialized. 

Correspondingly, the final spectra and elliptic flow 
show very little sensitivity to the initialization of 7r mn , 
as seen in Fig. 1181 With vanishing initial viscous pres- 
sure, viscous effects on the final flow anisotropy are a 
little weaker (dotted lines in Fig. [T8|) . but this difference 
is overcompensated in the total elliptic flow by slightly 
stronger anisotropies of the local rest frame momentum 
distributions at freeze-out (dashed lines in Fig. [TH]). For 
shorter kinetic relaxation times tv, the differences result- 
ing from different initializations of 7r mn would be smaller 
still. 



B. Kinetic relaxation time tv 

While the finite relaxation time for the viscous pres- 
sure tensor in the Israel-Stewart formalism eliminates 
problems with superluminal signal propagation in the 
relativistic Navier-Stokes theory, it also keeps the vis- 
cous pressure from ever fully approaching its Navier- 
Stokes limit ir mn = 2r\o mn . In this subsection we ex- 
plore how far, on average, the viscous pressure evolved 
by VISH2+1 deviates from its Navier-Stokes limit, and 
how this changes if we reduce the relaxation time tv by 
a factor 2. 

In Figure [TH] we compare, for central Cu+Cu colli- 
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FIG. 18: (Color online) Differential elliptic flow V2(pr) for 
pions from fe = 7fm Cu+Cu collisions with SM-EOS Q. Thick 
lines reproduce the pion curves from Figure [T2l obtained with 
7T mn = 2r]a mn at initial time to, while thin lines of the same 
type show the corresponding results obtained by setting ini- 
tially 7r mn = 0. 

sions, the time evolution of the scaled viscous pressure 
tensor, averaged in the transverse plane over the thermal- 
ized region inside the freeeze-out surface, with its Navier- 
Stokes limit, for two values of tv, ~ 3r/ / sT ~ r^ lass / 2 
and tv =T^ lass /4. For the larger relaxation time, the de- 
viations from the Navier-Stokes limit reach 25-30% at 
early times, but this fraction gradually decreases at later 
times. For the twice shorter relaxation time, the frac- 
tional deviation from Navier-Stokes decreases by some- 
what more than a factor 2 and never exceeds a value of 
about 10%. 

Figure I2TJ1 shows that, small as they may appear, these 
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FIG. 19: (Color online) Time evolution of the two in- 
dependent viscous pressure tensor components ■n TT and 
E — ■K xx +-K yy for central Cu+Cu collisions (solid lines), com- 
pared with their Navier-Stokes limits 2r\a TT and 2r](a xx +a yv ) 
(dashed lines), for two values of the relaxation time, 
T w = 3rj/sT (thick lines) and tv = 1.5n/sT (thin lines). All 
quantities are scaled by the thermal equilibrium enthalpy e+p 
and transversally averaged over the thermalized region inside 
the decoupling surface. 



deviations of 7r mn from its Navier-Stokes limit 2r)(j mn (es- 
pecially on the part of the decoupling surface correspond- 
ing to early times t— to) still play an important role for 
the viscous reduction of elliptic flow observed in our cal- 
culations. While a decrease of the relaxation time by a 
factor 2 leads to only a small reduction of the viscous sup- 
pression of flow anisotropies (dotted lines in Fig. [20]), the 
contribution to i^Cpt) resulting from the viscous correc- 
tion ~ VmVn^ mn to the final particle spectra is reduced 
by about a factor 2, too, leading to a significant over- 
all increase of V2{pr) in the region pT>lGeV/c. To 
avoid strong sensitivity to the presently unknown value 
of the relaxation time tv in the QGP, future extractions 
of the specific shear viscosity rj/s from a comparison be- 
tween experimental data and viscous hydrodynamic sim- 
ulations should therefore be performed at low transverse 
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FIG. 20: (Color online) Differential elliptic flow V2(pt) for n~ 
from b = 7fm Cu+Cu collisions with SM-EOS Q, calculated 
from viscous hydrodynamics with two different values for the 
relaxation time tv. Thick lines reproduce the pion curves 
from Figure 1121 thin lines show results obtained with a twice 
shorter relaxation time. For the standard (twice larger) classi- 
cal relaxation time value tv — 6rj/sT [2?], [HI deviations from 
ideal hydrodynamics would exceed those seen in the thick 
lines. 



momenta, pt < 1 GeV/c, where our results appear to be 
reasonably robust against variations of tv . 



C. Breakdown of viscous hydrodynamics at high px 

As indicated by the horizontal dashed lines in Figs. [5] 
and 1151 the assumption |<S/| <C |/ eq | under which the vis- 
cous hydrodynamic framework is valid breaks down at 
sufficiently large transverse momenta. For a quantita- 
tive assessment we assume that viscous hydrodynamic 
predictions become unreliable when the viscous correc- 
tions to the particle spectra exceed 50%. Fig. [8] shows 
that the characteristic transverse momentum p* T where 
this occurs depends on the particle species and increases 
with particle mass. To be specific, we here consider p^, 
for pions — the values for protons would be about 15% 
higher. The discussion in the preceding subsection of the 
Tv-dependence of viscous corrections to the final spec- 
tra makes it clear that reducing tv will also push p^ to 
larger values. Since we do not know tv we refrain from a 
quantitative estimate of this effect. 

In Fig. [2T] we show the breakdown momentum p* T for 
pions as a function of the peak initial energy density in 
the fireball center (i.e. indirectly as a function of colli- 
sion energy), for both central and peripheral Cu+Cu col- 
lisions. (The initial time was held fixed at ro = 0.6fm.) 
Generically, p^ rises with collision energy. The anomaly 
at low values of e(r=0) results, as far as we could as- 
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FIG. 21: (Color online) Characteristic transverse momentum 
where the viscous corrections to the final pion spectrum 
become so large (> 50%) that the spectrum becomes unreli- 
able, as a function of the initial energy density in the center 
of the fireball. Stars are for central Cu+Cu collisions, open 
circles for peripheral Cu+Cu collisions at b — 7 fm. Note that 
identical e(r=0) values correspond to higher collision energies 
in peripheral than in central collisions, values are higher 
for more massive hadrons (see Fig. [SJ , and they also increase 
for smaller relaxation times tv (see discussion of Fig. I20|l . 



certain, from effects connected with the phase transition 
in SM-EOS Q. The rise of p* T with increasing e(r=0) re- 
flects the growing fireball lifetime which leads to smaller 
viscous pressure components at freeze-out. This lifetime 
effect is obviously stronger for central than for peripheral 
collisions, leading to the faster rise of the stars than the 
open circles in Fig.[HU Taking further into account that a 
given beam energy leads to higher e(r=0) values in cen- 
tral than in peripheral collisions such that, for a given 
experiment, the peripheral collision points are located 
farther to the left in the figure than the central collision 
points, we conclude that in central collisions the validity 
of viscous hydrodynamics extends to significantly larger 
values of px than in peripheral collisions: Viscous effects 
are more serious in peripheral than in central collisions. 



VI. SUMMARY AND CONCLUSIONS 

In this paper, we numerically studied the shear vis- 
cous effects to the hydrodynamic evolution, final hadron 
spectra, and elliptic flow v 2 , using a (2+l)-dimensional 
causal viscous hydrodynamic code, VISH2+1, based on 
the 2nd order Israel-Stewart formalism. Using a fixed set 
of initial and final conditions, we explored the effects of 
shear viscosity for a "minimally" [19( viscous fluid with 
2 = in central and peripheral Cu+Cu collisions, com- 
paring the evolution with two different equations of state, 
an ideal massless parton gas (EOS I) and an EOS with a 
semirealistic parametrization of the quark-hadron phase 



transition (SM-EOS Q). Final hadron spectra and their 
elliptic flow were calculated from the hydrodynamic out- 
put using the Cooper-Frye prescription. 

We found that shear viscosity decelerates longitudi- 
nal expansion, but accelerates the build-up of transverse 
flow. This slows the cooling process initially, leading to 
a longer lifetime for the QGP phase, but causes acceler- 
ated cooling at later stages by faster transverse expan- 
sion. Viscous pressure gradients during the mixed phase 
increase the acceleration during this stage and slightly re- 
duce its lifetime. They counteract large gradients of the 
radial velocity profile that appear in ideal fluid dynamics 
as a result of the softness of the EOS in the mixed phase, 
thereby de facto smoothing the assumed first-order phase 
transition into a rapid cross-over transition. In the end 
the larger radial flow developing in viscous hydrodynam- 
ics leads to flatter transverse momentum spectra of the 
finally emitted particles, while their azimuthal anisotropy 
in non-central heavy-ion collisions is found to be strongly 
reduced. 

Although the viscous hardening of the hadron px- 
spectra can be largely absorbed by retuning the initial 
conditions, starting the transverse expansion later and 
with lower initial entropy density [3ll . |34| , this only acer- 
bates the viscous effects on the elliptic flow v 2 which in 
this case is further reduced by the decreased fireball life- 
time. The reduction of the elliptic flow v 2 by shear vis- 
cous effects is therefore a sensitive and robust diagnostic 
tool for shear viscosity in the fluid [62j . 

Our results indicate that in semiperipheral Cu+Cu col- 
lisions even a "minimal" amount of shear viscosity [l9[ 
causes a reduction of v 2 by almost 50% relative to ideal 
fluid dynamical simulations. In the present paper we ex- 
plored the origin of this reduction in great detail. The 
effects observed by us for Cu+Cu collisio ns 1371 are larger 
than those recently reported in Refs. [36l.l38l] for Au+Au 
collisions. While some of these differences can be at- 
tributed to an increased importance of viscous effects 
in smaller systems [58], the bulk of the difference ap- 
pears to arise from the fact that the different groups 
solve somewhat different sets of viscous hydrodynamic 
equations [H, [63|. (See also the recent interesting sug- 
gestion by Pratt [64| for a phenomenological modification 
of the Isreal-Stewart equations for systems with large ve- 
locity gradients.) This raises serious questions: if the- 
oretical ambiguities in the derivation of the viscous hy- 
drodynamic equations reflect themselves in large varia- 
tions of the predicted elliptic flow, any value of the QGP 
shear viscosity extracted from relativistic heavy-ion data 
will strongly depend on the specific hydrodynamic model 
used in the comparison. A reliable quantitative extrac- 
tion of rj/s from experimental data will thus only be pos- 
sible if these ambiguities can be resolved. 

Our studies show that shear viscous effects are 
strongest during the early stage of the expansion phase 
when the longitudinal expansion rate is largest. At later 
times the viscous corrections become small, although not 
negligible. Small non-zero viscous pressure components 
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along the hadronic decoupling surface have significant 
effects on the final hadron spectra that grow with trans- 
verse momentum and thus limit the applicability of the 
viscous hydrodynamic calculation to transverse momenta 
below 2-3GeV/c. depending on impact parameter, colli- 
sion energy and particle mass. Viscous effects are more 
important in peripheral than in central collisions, and 
larger for light than for heavy particles. They increase 
with the kinetic relaxation time for the viscous pressure 
tensor. Since the breakdown of viscous hydrodynamics is 
signalled by the theory itself, through the relative mag- 
nitude of the viscous pressure, the applicability of the 
theory can be checked quantitatively case by case and 
during each stage of the expansion. 

For the kinetic relaxation times tv considered in the 
present work, sensitivities to the initial value of the vis- 
cous pressure tensor were found to be small and practi- 
cally negligible. Sensitivity to the value of t„ was found 
for the hadron spectra, especially the elliptic flow, at 
large transverse momenta. This leads us to suggest to 
restrict any comparison between theory and experiment 
with the goal of extracting the shear viscosity rj/s to the 
region px < 1 GeV/c where the sensitivity to r n is suffi- 
ciently weak. 

The dynamical analysis of shear viscous effects on the 
momentum anisotropy and elliptic flow in non-central 
collisions reveals an interesting feature: The total mo- 
mentum anisotropy receives two types of contributions, 
the first resulting from the anisotropy of the collective 
flow pattern and the second arising from a local momen- 
tum anisotropy of the phase-space distribution function 
in the local fluid rest frame, reflecting viscous corrections 
to its local thermal equilibrium form. During the early 
expansion stage the latter effect (i.e. the fact that large 
viscous pressure effects generate momentum anisotropics 
in the local fluid rest frame) dominate the viscous effects 
on elliptic flow. At later times, these local momentum 
anisotropies get transferred to the collective flow profile, 
manifesting themselves as a viscous reduction of the col- 
lective flow anisotropy. The time scale for transferring 
the viscous correction to «2 from the local rest frame 
momentum distribution to the collective flow pattern ap- 
pears to be of the same order as that for the evolution of 
the total momentum anisotropy itself. 

Several additional steps are necessary before the work 
presented here can be used as a basis for a quantita- 
tive interpretation of relativistic heavy- ion data. First, 
the abovementioned ambiguity of the detailed form of 
the kinetic evolution equations for the viscous pressure 
must be resolved. Second, the equation of state must be 
fine-tuned to lattice QCD data and other available infor- 
mation, to make it as realistic as presently possible. The 
hydrodynamic scaling of the final elliptic flow vi with 
the initial source eccentricity e x [65j and its possible vio- 
lation by viscous effects need to be explored [43j | , in order 
to assess the sensitivity of the scaled elliptic flow V2/e x 
to details of the model used for initializing the hydrody- 
namic evolution [l4[ ■ The temperature dependence of the 



specific shear viscosity Tjls, especially across the quark- 
hadron phase transition [60l . l6l| , must be taken into ac- 
count, and bulk viscous effects, again particularly near 
T c , must be included. To properly account for the highly 
viscous nature of the hadron resonance gas during the 
last collision stage it may be necessary to match the vis- 
cous hydrodynamic formalism to a microscopic hadronic 
cascade to describe the last part of the expansion until 
hadronic decoupling [60| . We expect to report soon on 
progress along some of these fronts. 

Note added: Just before submitting this work for pub- 
lication we became aware of Ref. [68[ where the form 
of the kinetic evolution equations for the viscous pres- 
sure is revisited and it is argued that Eqs. (|4l5p must 
be amended by additional terms which reduce the strong 
viscous suppression of the elliptic flow observed by us 
[(33| . While details of the numerical results will obviously 
change if these terms are included (c/. Rcfs. [H, HH), 
our discussion of the driving forces behind the finally ob- 
served viscous corrections to ideal fluid results and of the 
evolution of these corrections with time is generic, and 
the insights gained in the present study are expected to 
hold, at least qualitatively, also for future improved ver- 
sions of VISH2+1 that prop erly take into account the 
new findings reported in [68| |. 
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APPENDIX A: EXPRESSIONS FOR n mn AND a mn 



The expressions for jt mn and a mn in Eq. (|T0|) are 
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APPENDIX B: VELOCITY FINDING 



As shown in [28[ , since we evolve all three components 
7r TT , tt tx , and 7r ry (one of which is redundant due to 
the constraint ir Tm u m = 0), the flow velocity and energy 
density can be found from the energy-momentum tensor 
components with the same efficient one-dimensional zero- 
search algorithm employed in ideal hydrodynamics [67l ]. 
This is important since this step has to be performed 
after each time step at all spatial grid points in order to 
evaluate the EOS p(e). 

Using the output from the numerical transport al- 
gorithm, one defines the two-dimensional vector M = 
(M x ,M y ) = (T TX -n TX , T T y~n T 'y). This is the ideal fluid 
part of the transverse momentum density vector; as such 
it is parallel to the tranverse flow velocity v± = (v x , v y ) 
Introducing further Mq = T rT — 
ergy density as 



one can write the en- 
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M = M Q - v±M, 
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where uj_ = 
M = 



-Vy is the transverse flow speed and 



M 2 +M 2 . One sees that solving for e requires 

only the magnitude of v \ which is obtained by solving 
the implicit relation [28|, [67j 
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by a one-dimensional zero-search. The flow velocity com- 
ponents are then reconstructed using 
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Note that this requires direct numerical propagation of 
all three components (ir TT , tt tx and ir Ty ) since the flow 
velocity is not known until after the velocity finding step 



has been completed. Hence the transversality constraint 
■n Tm u m = cannot be used to determine, say, ■k tt from 
tt tx and tt tv . However, it can be used after the fact to 
test the numerical accuracy of the transport code. 



APPENDIX C: n mn IN TRANSVERSE POLAR 
COORDINATES 

Although VISH2+1 uses Cartesian (x,y) coordinates 
in the transverse plane, polar (r, <p) coordinates may 
be convenient to understand some of the results in the 
limit of zero impact parameter where azimuthal sym- 
metry is restored. In (t, r, <p, rj) coordinates the flow 
velocity takes the form u m = Y_i_(l, v r , v^, 0), with 

7_i_ = 1/ y/T— ~v\ = 1/ ^Jl—v 2 —r 2 v 2 . The polar coordinate 

components of the shear pressure tensor components ir mn 
are obtained from those in (r, x, y, rj) coordinates by the 
transformations 
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with cos 4> = x/r and sin<f) ~y/r. In terms of these the 
independent components S and A of Eqs. (|19I20[) are 
given as 

S = 7r rr + r 2 7r^, 

A = cos(2(/))(7r rr -r 2 7r^) - 2sin(2(/.)r7r r *, (C2) 

from which we easily get 

2tt xx = 7r rr ((l + cos(20)) +r 2 7r^((l-cos(20)), 
2lr yy = ^rr^ _ cos ( 2 0)) + r 2 7T^((l + cos(20))[C3) 

Note that azimuthal symmetry at b = implies ir r ^ = 
and a vanishing azimuthal average for A: (A)^ = or 



APPENDIX D: TESTS OF THE VISCOUS 
HYDRO CODE VISH2+1 

1. Testing the ideal hydro part of VISH2+1 

When one sets ir mn = initially and takes the limit 
?/ = 0, VISH2+1 simulates the evolution of an ideal fluid, 
and its results should agree with those of the well-tested 
and publicly available (2+l)-dimensional ideal fluid code 
AZHYDRO [4l|. Since VISH2+1 was written indepen- 
dently, using only the flux-corrected SHASTA transport 
algorithm from the AZHYDRO package [U [66| in its 
evolution part, this is a useful test of the code. The 
left panel in Fig. |2"21 shows that, for identical initial and 
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FIG. 22: (Color online) Left: Differential elliptic flow V2(pr) for n~ from 6 = 4fm Cu+Cu collisions and b — 7im Au+Au 
collisions, using EOS Q. Results from VISH2+1 for rj = and n mn = (dashed lines) are compared with the ideal fluid code 
AZHYDRO (solid lines) . Right: V2 (pr) for -k~ from Cu+Cu collisions at impact parameters 6 = 4 and 7 fm, comparing VISH2+1 
evolution with EOS Q (dashed) and SM-EOS Q (solid) in the ideal fluid limit r\ = 0, 7r mn = 0. 



final conditions as described in Sec. HH the two codes 
indeed produce almost identical results. The small dif- 
ference in the Au+Au system at b = 7 fm is likely due 
to the slightly better accuracy of AZHYDRO which, in 
contrast to VISH2+1, invokes an additional timesplitting 
step in its evolution algorithm. 

When comparing our VISH2+1 results with AZHY- 
DRO we initially found somewhat larger discrepancies 
which, however, could be traced back to different versions 
of the EOS used in the codes: EOS Q in AZHYDRO, 
the smoothed version SM-EOS Q in VISH2+1. In the 
left panel of Fig. [22] this difference has been removed, by 
running also VISH2+1 with EOS Q. In the right panel we 
compare VISH2+1 results for EOS Q and for SM-EOS Q, 
showing that even the tiny rounding effects resulting from 
the smoothing procedure used in SM-EOS Q (which ren- 
ders the EOS slightly stiffer in the mixed phase) lead to 
differences in the elliptic flow for peripheral collisions of 
small nuclei which exceed the numerical error of the code. 
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FIG. 23: (Color online) Comparison between the analyti- 
cal temperature evolution for (0+I)-d boost-invariant Navier- 
Stokes viscous hydrodynamics (solid line) and numerical re- 
sults from VISH2+I with homogeneous transverse initial en- 
ergy density profiles (dashed line). 



2. Comparison with analytical results for (0+l)-d 
boost-invariant viscous hydro 



ing analytic solution for the temperature evolution [11 



For boost-invariant longitudinal expansion without 
transverse flow, the relativistic Navier-Stokes equations 
read Fill 
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For an ideal gas EOS p=\e^T A this leads to the follow- 
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To test our code against this analytical result we initialize 
VISH2+1 with homogeneous tranverse density distribu- 
tions (not transverse pressure gradients and flow) and 
use the Navier-Stokes identification ir rnn — 2r\a mn in the 
hydrodynamic part of the evolution algorithm, sidestep- 
ping the part of the code that evolves 7r" m kinetically. 
It turns out that in this case the relativistic Navier- 
Stokes evolution is numerically stable. Fig. |2"31 compares 
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the numerically computed temperature evolution from 
VISH2+1 with the analytic formula (|D3|) . for ry/s^ 0.08 
and Tq = 360 McV at r = 0.6 fm/c. They agree perfectly. 



3. Reduction of VISH2+1 to relativistic 
Navier-Stokes theory for small r\ and tv 

Having tested the hydrodynamic part of the evolution 
algorithm in Appendix ID 11 we would like to demon- 
strate also the accuracy of the kinetic evolution algorithm 
that evolves the viscous pressure tensor components. A 
straightforward approach would be to take VISH2+1, set 
the relaxation time tv as close to zero as possible, and 
compare the result with a similar calculation as in Ap- 
pendix Q7T] where we sidestep the kinetic evolution algo- 
rithm and instead insert into the hydrodynamic evolution 
code directly the Navier-Stokes identity n mn = 2 r qa mn . 
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FIG. 24: (Color online) Differential elliptic flow V2(pr) for 
gluons from b = 7 fm Cu+Cu collisions, calculated with ideal 
hydrodynamics (blue dashed line), relativistic Navier-Stokes 
(NS) hydrodynamics (light blue lines), and Israel-Stewart 
(IS) viscous hydrodynamics with ^ = 2 ^ v and tw = 0.03 fm/c 
(red lines), using EOS I. The lines for NS and IS viscous hy- 
drodynamics are almost indistinguishable. Solid lines show 
the full results from viscous hydrodynamics, dotted lines ne- 
glect viscous corrections to the spectra and take only the flow 
anisotropy effect into account. 

Unfortunately, this naive procedure exposes us to the 
well-known instability and acausality problems of the rel- 
ativistic Navier-Stokes equations. The suggested proce- 
dure only works if a set of initial conditions and transport 
coefficients can be found where these instabilities don't 
kick in before the freeze-out surface has been reached. 

We found that sufficiently stable evolution of the rel- 
ativistic Navier-Stokes algorithm (i.e. of VISH2+1 with 
the identification ir mn = 2r\a mn ) can be achieved for stan- 
dard initial density profiles in Cu+Cu collisions and the 
simple ideal gas equation of state EOS 1 by choosing 
a very small and temperature dependent specific shear 
viscosity ^ = 0.01 2 oomcV = 2 Gey ^ or ^ e I sraei_ Stewart 



evolution we use a relaxation time which is correspond- 
ingly short: tv = = 0.03 fm/c. 

Figure [M] shows the differential elliptic flow V2(pt) for 
gluons in 6 = 7 fm Cu+Cu collisions evolved with these 
parameters. The dashed line gives the ideal fluid result. 
The solid and dotted lines show the total elliptic flow and 
the anisotropic flow contribution to ^(pt), respectively, 
similar to the left panel Fig. [T3J There are two solid 
and dotted lines with different colors, corresponding to 
Israel- Stewart and Navier-Stokes evolution; they are in- 
distinguishable, but clearly different from the ideal fluid 
result. We conclude that, for small shear viscosity rj/s 
and in the limit tv — > 0, the second-order Israel-Stewart 
algorithm reproduces the Navier-Stokes limit and that, 
therefore, VISH2+1 evolves the kinetic equations for ir mn 
accurately. 



APPENDIX E: HYDRODYNAMICS VS. BLAST 
WAVE MODEL 

As discussed in Sec. lIIIBl the viscous corrections to the 
final pion spectra from the hydrodynamic model have a 
different sign (at least in the region Vj_> 1 GeV) than 
those originally obtained by Teaney In this Ap- 

pendix we try to explore the origins of this discrepancy. 
We will see that the sign and magnitude of viscous correc- 
tions to the (azimuthally averaged) particle spectra are 
fragile and depend on details of the dynamical evolution 
and hydrodynamic properties on the freeze-out surface. 
Fortunately, they same caveat does not seem to apply 
to the viscous corrections to elliptic flow where hydro- 
dynamic and blast wave model calculations give qualita- 
tively similar answers. 

Following Teaney's procedure, we calculate ir mn in the 
Navier-Stokes limit -K mn = 2rja mn . We do this both in 
the blast wave model and using the results for a mn from 
VISH2+1. For the blast wave model we assume like 
Teaney freeze-out at constant r with a box-like density 
profile e(r) = ed ec 0(Ro—i"), where edec = 0.085 GeV/fm 3 
is the same freeze-out energy density as in the hydro- 
dynamic model for EOS I, and Ro = 6 fm. The veloc- 
ity profile in the blast wave model is taken to be lin- 
ear, u r (r) — a,Qj^9(R —r), with a Q = 0.5; freeze-out is 

assumed to occur at Td cc = 4.1 fm/c. Ro, ao and Td cc 
are somewhat smaller than in Ref. [l6| since we study 
Cu+Cu instead of Au+Au collisions. We concentrate 
here on a discussion of ir rr for illustration; the expression 
for <r rr is found in Ref. [H, Eq. (Allc). While ir rr from 
VISH2+1 differs from 2r\a rr due to the finite relaxation 
time tv (see Sec. IV Cp . we have checked that the signs 
of these two quantities are the same on the freeze-out 
surface so that our discussion provides at least a qualita- 
tively correct analysis of the viscous spectra corrections 
in the two models. 

In Fig. [25] we compare the freeze-out profiles for the ra- 
dial flow velocity and 2r]o rr from the blast wave model. 
In spite of qualitative similarity of the velocity profiles, 
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FIG. 25: (Color online) Top row: Velocity profiles from the blast wave model (left) and from the hydrodynamic model with EOS I 
at fixed times (middle) and along the decoupling surface (right). Bottom row: The corresponding profiles for the transverse 
shear viscous pressure 7r rr in the Navier-Stokes limit, 7r rr = u . Calculations are for central Cu+Cu collisions, and the 

curves in the middle panels correspond to the times r = l, 2, 4, and 6 fm/c. See text for discussion. 



the freeze-out profiles of 2r/a rr arc entirely different and 
even have the opposite sign in the region where most of 
the hydrodynamic particle production occurs (left and 
right columns in Fig. [25]). The middle column shows 
that at fixed times r, the hydrodynamic profile for 2r/a rr 
shows some similarity with the blast wave model in that 
2r\a Tr is positive throughout most of the interior of the 
fireball. What matters for the calculation of the spec- 



tra via Eq. (|12p . however, are the values of 2r\o rr on the 
freeze-out surface £ where they are negative, mostly due 
to radial velocity derivatives. This explains the opposite 
sign of the viscous correction to the spectra in the hydro- 
dynamic model and shows that, as far as an estimate of 
these viscous corrections goes, the blast wave model has 
serious limitations. 
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